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The Basel problem consists in evaluating the series 
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It was tackled by many mathematicians in the 17th century, most notably Pietro Men- 
goli, Gottfried Wilhelm Leibniz and Jakob Bernoulli jT, p42]. It was finally solved by 
Leonhard Euler in 1735. He showed that 
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a result that still astonishes everyone who sees it for the first time. 

Even though Euler’s first proof was not rigorous, it convinced contemporary mathe¬ 
maticians. Since then, many proofs of this identity have been found whose rigor stands 
up to the scrutiny of modern mathematics. A (non-exhaustive) list of such proofs can 
be found in [3]. Some of them are widely praised for their esthetic value, see [TJ Chapter 
8 ]- 

One proof relies on the series expansion 
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Setting h = 1 in this expansion, combined with some elementary but nontrivial alge¬ 
braic manipulations, yields the desired result. However, that approach only relocates 
the difficulty because proving this expansion is not an easy task. Instead of tackling 
this problem head on, we will turn our attention to a completely different branch of 
mathematics. 

Modified equations are a powerful tool from numerical analysis. In the next section, we 
will introduce this concept and start studying an example. As the discussion progresses, 
it will gradually become clear that the modified equation in this example is intimately 
related to the series expansion ©• 
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The main goal of this article is to show off the concept of modified equations, a beau¬ 
tiful idea which rarely travels beyond the borders of numerical mathematics. However, 
the point of this work is not to provide a general overview of this subject. Excellent 
review texts are already available, for example [5j Chapter IX]. The point of this article 
is to show that modified equations have potential use outside of numerical analysis. To 
illustrate this, at the end of this paper we will use a modified equation to prove Equation 
m and by corollary the Basel problem. 


Modified equations 

An important technique to study the behavior of a numerical method for an ODE is 
backward error analysis. Instead of comparing a discrete solution to a solution of the 
original ODE, one looks for a modified equation whose solutions exactly interpolate 
the numerical solutions. To analyze the behavior of the numerical integrator one then 
compares the modified equation with the original one. In general a modified equation 
involves a formal power series, so to extract information in a rigorous way one needs to 
truncate it. In many cases useful estimates for the truncation error are available, and 
very strong results can be proved using modified equations. The most famous one is 
the fact that symplectic integrators for Hamiltonian systems nearly conserve the energy 
over very long time intervals. 

Let us illustrate the concept of a modified equation in the context of the linear differ¬ 
ential equation 

x(t) = —x(t), 

where the dot denotes differentiation with respect to t. The dependence of x on t will 
often be suppressed in what follows. We consider the explicit Euler discretization with 
step size h of this differential equation, which is obtained by replacing the derivative of 
x by the finite difference dj+1 h — and x itself by Xj. For each integer j the quantity Xj 
is intended to be an approximation of x{jh). This gives us the difference equation 

Xj -|_i — Xj = —hxj. 

Its modified equation is a differential equation x = Fh(x ) whose solutions satisfy the 
difference equation, in the sense that x(t + h) — x(t) = —hx(t). The right hand side of 
the modified equation is a function F : R + xR->K: (h, x) H > Fh(x), which we assume 
to have a power series expansion 

Fh(x) = f 0 (x) + hfi(x) + h 2 f 2 (x ) + .... 

We will assume that this series converges for the step size h under consideration. In 
most cases the power series defining a modified equation does not converge for any h, 
but in this paper we will only encounter examples where convergence takes place. 

The coefficients fj(x) of the power series can be determined as follows. Using Taylor 
expansion we can write the condition that x(t + h) — x(t) = — hx(t ) as 

hx(t) -\ —- x(t ) H- x^ 3 \t) + ... = —hx(t). 
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The higher derivatives can be written as 
x = F} l (x)F h (x) 

= ( fo(x ) + hf[(x) + h 2 f 2 {x) + ...) (f 0 (x) + hfi(x) + h 2 f 2 {x) + ...), 
x (3) = Ffi(x)F h (x) 2 + F' h {x) 2 F h {x ) 

= (/o0*0 + hf"(x) + h 2 f"(x) + ...) (/ 0 (x) + hfi(x) + h 2 f 2 (x ) + .. .) 2 

+ (/o(z) + hf[(x) + h 2 f' 2 {x) + .. .) 2 (f 0 (x) + hfi(x) + h 2 f 2 (x) + ...), 

and so on, where the prime denotes differentiation with respect to x. Hence the condition 
becomes 

h (/o + hh + h 2 f 2 + ...) 

+ ~ (/o/o + hf[fo + hfoh + ■■■) + — (/o/o + fofo + •••) + ••• = — hx, 
where the argument x of the fi has been omitted. Grouping terms by order in h we find 

Hfo + x) + h~ ^/i + -fofoj + h 3 ^/ 2 + -f[fo + ^/o/i + g/ 0/0 + g/o 2 /o^ + ■ • • = 0. 

For a power series to be equal to zero, all of the coefficients must be zero. From the first 
order term we find that fo(x) = —x. From the second order term we then learn that 
f\ (x) = —\x, and from the third order term that f 2 (x) = —\x. Hence the modified 
equation is 

h h 2 

x = —x - x - x — . . . . 

2 3 

We see that the leading order term agrees with the original differential equation. The 
higher order terms reflect the discretization error. 

Since the difference equation Xj +1 = xj — hxj is linear it can be solved exactly and the 
modified equation doesn’t provide any new information. However, the above procedure 
can be applied just as well to nonlinear difference equations. The only price we pay is 
that in the nonlinear case the terms tend to become more and more complicated as the 
order in h increases. 

For first order ODEs the notion of a modified equation is well-established [5j Chapter 
IX]. It is easily extended to second order equations [9], which is the relevant setting for 
this paper. 

Definition. Let / : M + x M 2 —)• R : (h,x,v) i->- Fh(x,v) and T : R + x R 3 —)• R : 
( h,x,y,z ) T/i( x,y,z) be smooth functions. The differential equation x = Ff l (x,x) 
is a modified equation for the second order difference equation ^>h(xj-i, Xj, ®j+i) = 0 
if for all sufficiently small h > 0, every solution x of x = E) l (x,x) satisfies SPh(x(t — 
h),x(t),x(t + h)) = 0 for all f £ R. 

In other words, for every solution x of the modified equation, the discrete curve 
(x(to + jh))j£z solves the difference equation. Usually the modified equation is ob¬ 
tained as a formal power series in the step size h that does not converge. To be rigorous 
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in general, the definition should be adapted to handle a right hand side Ff l (x, x) that is 
a formal power series, but in this paper we will only encounter convergent series. 

As an example, consider the differential equation x = — x and its Stormer-Verlet 
discretization, which is obtained by replacing the second derivative of x by the finite dif¬ 
ference Xj+1 m This is a famous example of a geometric (i.e. structure-preserving) 

numerical integrator [B|. It gives us the difference equation 

Xj .|_i — 2 Xj + Xj -1 = —h 2 Xj. (2) 

In this particular example the modified equation is unusually simple. It turns out to be 
sufficient to look for a modified equation of the form 

x = F h {x) = f 0 (x) + h 2 f 2 (x) + h 4 f 4 (x) + .... (3) 


In general we should also include odd order terms and let the /j depend on x as well. It 
follows from Equation ([3]) that 

x (3) = F^x, 

= FHx 2 + FiF h , 
x {5) = F^x 3 + 3 Fl[F h x + Ftfx, 

x (6) = if } x 4 + 6 F^F h x 2 + 5 F^Fix 2 + 3 F^F 2 + F’ 2 F h , 


where the argument x of Fh has been omitted. We identify x(t) = Xj and Xj±\ = x(t±h). 
Using Taylor expansion we find 


. . /i 2 h 3 

a: o+i = x ± hx 4- x ± —x 

3 2 6 
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*L X («) ± 2L, 
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.(7) 


+ .... 


Plugging this into the difference equation © and replacing derivatives using ([2]) and fU 
we find 


— h 2 x = h 2 x + —x^ + x^ + ... 

12 360 

= h 2 (f 0 + h 2 h + h A U) + ^ (/"x 2 + h 2 f%x 2 + /'/ 0 + h 2 .flf2 + h 2 f'f 0 ) 

+ (/o + 6/o ^ fox 2 + 5/o fox 2 + 3/q/o + fo /o) + ■■■ 

= h 2 f 0 + h A (f 2 + ^ (/"x 2 + 

+ + Yo /0/2 + /2/0) 

+ 360 + 6 ^d^/cU 2 + 5/o/ qX 2 + 3/q/o + /o 2 /o) ^ + — 
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Figure 1: The harmonic oscillator. Plotted are an exact solution (dashed curve), a 
solution of the Stormer-Verlet discretization with step size h = 1 (dots), a solution of 
the modified equation up to order two (light solid curve), and a solution of the modified 
equation up to order four (dark solid curve). 


The /) 2 -terrn of this equation gives us fo(x) = —x, hence f( = —1 and f(( = 0. The 
h 4 - term then reduces to f 2 (x) = — hence f 2 = — py and f'f = 0. Finally, the /i 6 -term 
gives 

x 

90 


, , s 1 ( X x \ 

/4 ( X) + 


360 


Therefore, the modified equation is 


h' 2 


h 4 


x = —x - x - x — ... 

12 90 


( 5 ) 


Note that the coefficients fj are uniquely determined by the above procedure. This 
shows that if there exists a modified equation that can be represented as a power series 

x = fo(x,x) + hfi(x,x) + h 2 f 2 (x,x) + h 4 f 3 (x,x) + ... 


then it is unique. 

Figure Q] shows a solution of the modified equation truncated after the /i 4 -term. We see 
that it agrees very well with the discrete flow, despite the large step size. Furthermore, 
we observe that the Stormer-Verlet discretization of the harmonic oscillator with step 
size h = 1 is 6-periodic. Since the difference equation has solutions with a period of 
exactly 6, so does the modified equation. This suggests that the modified equation is 
x = —^x. 

Because the difference equation d2J) is linear, this observation is easy to prove. Indeed, 
we can find an explicit expression for the general solution of the difference equation: 

Xj = Ae~ 2ije + Be 2ije , 

where 9 = arcsin(^). Hence we can always find an interpolating curve of the form 
x(t) = Ae~ 2lte / h + Be 2lte / h . This curve satisfies the differential equation 



2 

If we set h = 1 this becomes x = —^x. 

At this point the connection between modified equations and Equation 0 begins to 
reveal itself. Our proof of Equation JT]) - and by corollary of the Basel problem - will 
be based on a comparison of the expressions © and © of the modified equation. The 
main difficulty is to find a general expression for the coefficients of the power series in 
Equation ©. 


5 




Determining the coefficients of (ED 


Another consequence of the linearity of the difference equation (J2]) is that the sum 
Xj_k + Xj + k can be written as Xj multiplied by a polynomial of degree 2k in h. More 
precisely, any solution of d2J) satisfies 


Xj—k T Xj-\-k 


2 T k 





where T k denotes the fc-th Chebyshev polynomial of the first kind (proof in the Ap¬ 
pendix). This implies that any solution of the modified equation satisfies 

/ h 2 

x(t — kh) + x(t + kh ) = 2T k ( 1 —— 

= (— l) k h 2k x(t) + terms of lower order in h. 



Now we are in a position to derive an explicit expression for the coefficients of the 
modified equation @0 Fix an arbitrary smooth curve x. For every j and k there holds 


x(t-jh)-2x(t)+x(t+jh) = (, jh) 1 2 x(t) + x (4) (f) + .. ■+ 2 |^| ! x^ 2k) (t) + 0{h 2k+2 ), 
or, in matrix form, 


/ x(t — h) — 2 x(t) + x(t + h) \ 


(1 1 ... 1\ 


( h 2 x{t) ^ 

x(t — 2 h) — 2 x(t) + x(t + 2 h) 

— 

2 2 2 4 2 2k 



\x(t — kh) — 2 x(t) + x(t + kh) / 


\k 2 A: 4 ... k 2k j 




+ 0{h 2k+2 ). 


Using Cramer’s rule we solve the above system of linear equations for h?x{t ), 


h 2 x(t) 


x{t — h) — 2 x(t) + x(t + h) 

1 . 

. 1 

x(t — 2 h) — 2 x(t) + x(t + 2 h) 

2 4 . 

. 2 2k 

x(t — kh) — 2 x(t) + x(t + kh) 

/c 4 . 

. k 2k 


1 1 ... 1 

2 2 2 4 2 2k 




k 2 fc 4 ... k 2k 




+ 0{h 2k+2 ). 


In the denominator we have a Vandermonde determinant that equals 

(to 2 n o' 2 -' 2 )- 

l<i<j<k 

1 The argument here is inspired by the derivation in [7] of suitable coefficients for finite difference 

methods. 
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The determinant in the numerator is more difficult to evaluate, but if we restrict our 
attention to the /i 2fc -term it becomes just as easy. Assuming the curve x is a solution of 
the modified equation, we can use Equation 0 to find that this term is 


x(t ) 


0 

1 

1 


1 

2 4 

1 

2 2k 

0 

2 4 

2 2k 


0 

(-1 ) k x(t) 

(k- l) 4 • 
k A 

. (k - l) 2fc 
k 2k 


(*- l) 4 • 

. (k - l) 2k 


-(k— i)! 4 1 n (/ 


— I 


x(t). 


x(t), 


Hence we find that the h 2k - term of h?x{t) equals 

which simplifies to 

k(2k — 1)! 1 J ' 

This is the desired explicit expression for the h 2k ~ 2 -term of Equation Q. Hence the 
modified equation of the difference equation xj + \ — 2Xj + Xj- \ = —h 2 Xj is 


x = — 


f ] 2{k Z 1 }' 2 h 2k ~ 2 x. 


k =1 


(2 k)\ 


( 8 ) 


Fitting the pieces together 

Equations d6j) and (JHJ) provide two expressions for the modified equation of the difference 
equation ©■ Since the modified equation, written in the form “x = power series” is 
unique, it follows that both expressions coincide: 
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2 k=v 


This expansion is relatively well-known. However, the proofs usually found in the liter- 
aturfU are quite complicated. Plugging in h = 1 we find 
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( 9 ) 


2 e.g. the ones sketched in [2J p. 384-385] and [8], p. 271]. 


7 















By elementary but nontrivial calculations (presented in the Appendix, following [8J p. 
265-266]) one can show that 


V 1 

^ k 2 ^ (2k)\ 

k =1 fc=l v ' 


which leads to the conclusion that 


oo .. o 

E l _ 7T 

P “ ~6~' 

k =1 

It is worth mentioning that the closed-form expression ([6]) of the modified equation 
in terms of the arcsine is not needed to arrive at this result. Indeed one can derive 
Equation Q from Equation ([8]) and the fact that solutions of the difference equation 
([2]) are 6-periodic for h = 1. The serendipitous observation of this periodicity was the 
inspiration for our unconventional solution of the Basel problem. 


Conclusion 

Modified equations are an important tool in numerical analysis. One can think of their 
study as reversing the discretization: one looks for a differential equation for which the 
discretization would have been exact. This continuous system then provides information 
about the discrete one. 

Usually, modified equations are given by formal power series. When for a specific 
example the power series does converge, the properties of the discrete system can be 
useful to evaluate it. This reverses the direction of thinking once more: we are back to 
using the discrete system to learn about the continuous one. However, in this context 
the discrete system is not a mere approximation of some continuous object, but leads us 
to an exact evaluation of a power series. 

At the moment, the scope of this method is limited to the example presented above. 
Nevertheless, it is very pleasing to find a connection between a problem that has fasci¬ 
nated mathematicians for centuries and the relatively new concept of modified equations 
for numerical integrators for ODEs. 
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Appendix 


Proposition 1. Let T k denote the fc-th Chebyshev polynomial of the first kind. Any 
solution of the difference equation Xj + \ — 2 Xj + xj- \ = —h 2 Xj satisfies 

h 2 \ 


Xj—k -|- Xj+k 2Tfc ^1 ^ J x j- 

Proof. As observed in the main text, solutions of the difference equation are of the form 

xj = Ae~ 2ij0 + Be 2ije , 

where 6 = arcsin ( f ). Therefore, 

Xj- k + x j+k = Ae~ 2i{j - k)e + Be 2i ^~ k)e + Ae~ 2i ^ +k)e + Be 2 ^ j+k)e 
= Ae~ 2i i e (e 2ifce + e~ 2ikd ) + Be 2i ^ e (e' 

= 2cos(2 k6)xj. 

Note that 


„-2 ik9 _|_ g 2 ik9 


cos 2(9 = 1 - 2 sin 2 0 = 1-2 


= 1 - 


h 2 


2 ) 2 ’ 

hence, by the trigonometric definition of the Chebyshev polynomials, 

h 2 


Proposition 2. 


cos(2 kO) = T k ( 1 - — 


OO 1 oo 

Et2= 3 E 
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k 2 (2 k)\ ' 
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Proof. We follow [8j p. 265-266]. Observe that for any k and l 
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k 2 k{k + 1) ife(fc + l)(fc + 2) k(k + l)...(k + l + l) 
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